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ABSTRACT 

We extend the one-dimensional, two-zone models of long-term protostellar disk evolution 
with infall of Zhu et al. to consider the potential effects of a finite viscosity in regions where the 
ionization is too low for the magnetorotational instability (MRI) to operate (the "dead zone"). 
We find that the presence of a small but finite dead zone viscosity, as suggested by simulations of 
stratified disks with MRI-active outer layers, can trigger inside-out bursts of accretion, starting 
at or near the inner edge of the disk, instead of the previously-found outside-in bursts with zero 
dead zone viscosity, which originate at a few AU in radius. These inside-out bursts of accretion 
bear a qualitative resemblance to the outburst behavior of one FU Ori object, V1515 Cyg, in 
contrast to the outside-in burst models which more closely resemble the accretion events in FU 
Ori and V1057 Cyg. Our results suggest that the type and frequency of outbursts are potentially 
a probe of transport efficiency in the dead zone. Simulations must treat the inner disk regions, 
R < 0.5 AU, to show the detailed time evolution of accretion outbursts in general and to observe 
the inside-out bursts in particular. 

Subject headings: accretion disks, stars: formation, stars: pre-main sequence 



1. Introduction 



The (re)discovery of the magnetorotational instability (MRI: e.g.. iBalbus & Hawleylll998l . and ref- 
erences therein), appears to resolve the long-standing problem of the ano malous viscosity in sufficie ntlv 
ionized accretion disks. To a crude approximation, this valid ates the use of IShakura & Sunvaevi (119731) "« 



viscosity" disks, though there are differences in detail (e.g.. iBalbus & Papaloizoulll999l : 



Gammie 



19961) . 



Constant a disks have been emp loyed in many situations, incl uding the evolution of pre-main sequence 
disks (e.g.. iHartmann et al.lll998h . However, as pointed out by iGammid (119961 ). thermal ionization levels 
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in protostellar and protoplanetary disks generally are so low that it is unUkely that the MRI operates ev- 
erywhere. Gammie suggested that transport in these regions might be limited to surface "active" layers, 
in which non-thermal ionization (cosmic rays, stellar X-rays) could allow the MRI to operate, while in the 



enabled beyond simple (quasi-steady) viscous disks ( 


Vorobvov & Basu 


2005, 


2006, 


2007 


2009: 


Zhu et al. 


2010a 


b; 


Martin & Lubow 


2011; 


Martin et al. 


2012a 1 


J). 



Gammid (119961) noted that if the MRI-active layer has a roughly constant surface density, the result 
is a pileup of material in the inner disk which eventually could beco me gravitationally unstable a nd result 



1996). While 



in a rapid burst of accretion, perhaps producing an FU Ori outburst (iHartmann & Kenvon 
early models of FU Qri accretion e vents relied on traditional thermal instability theory (IClarke et al.lll990 ; 
Bell&Linlll994h . IZhu et al.l (|2007n showed that the rapidly-accreting region of FU Ori extends much further 
in radius from the central star than can be achieved with this theory. The very large amount of mass accreted 
m one of FU Ori's outbursts (~ 10"^ Mq) over such short timescales (~ 10^ yr) requires a large amount 
of material to be present in the disk at relatively s mall radii. This is a n a tural result of models in which 

tri ggers the outbursts (Armitage et al. 200 ll: \^robyov & Basu 2006, 2007 . 



gravitational instability (GI) t ri ggers the outbursts (jArmitage et al 
2OO9I : Izhu et al.ll20093 . boiOal lbl- Matin & Lubowlbol ll : iMartinetal. 



2012al lbh. 



Zhu et al.l (|2010br) presented one-dimensional, two-layer evolutionary disk models including infall from 
a rotating protostellar cloud, and showed that they could qualitatively reproduce the main features of the 
outbursts of FU Ori and V1057 Cyg, with rapid rise times and slowly-decaying accretion. The simulations 
included irradiation from the central star, but did not take into account the accretion luminosity, which dur- 
ing outbursts can be far larger than the st ellar photospheric radiation. While the geometry of disk accretion 
may not favor self-irradiation ( Belllll999 ). during infall the dusty opaque envelope will ac t as a blanket, rera- 
diating a significant portion of the accretion luminosity toward the disk. In addition the IZhu et al.l (1201 Obi) 
calculations assumed that the central layer - the dead zone - had no viscosity unless it became gravitation- 
ally unstable. However, detailed shearing box simulations of the MRI in disks with a stratified structure 
and a resistivity that increases toward the midplane indicate that MHD turbulenc e generated in the upper 



MRI-active layers produces some hydrodynam ic turbulence in the inactive layers (|Fleming & Stonel 12003 



Okuzumi & Hirosell201ll : iGressel et al.ll2012|) . Furthermore, these investigations suggest that this turbu- 
lence creates Maxwell stresses, which result in a small but non-zero viscosity that can transport angular 
momentum outward and thus mass inward. 

These considerations motivate further investigations of protostellar accretion using the two-layer one- 
dimensional model. We find that the position and properties of the inner boundary, the inclusion of ir- 
radiation by the accretion luminosity generated in the inner disk, and any non-zero dead zone viscosity 
have significant effects on the resulting bursts of mass accretion. Our treatment is sufficiently limited to 
preclude detailed predictions, but the qualitative behavior is suggestive, given that continuing observations 
of protostars and pre-main sequence stars are increasingly found to exhibit a wide variety of accretion 



events, beyond the large FU Ori outbursts ( Hartmann & Kenyon 



Reipurth & Aspinll2010l : lAspin et aLll2010l : 



Covev et al. 



2011 



1996; 



Lorenzetti et al. 



Herbigl2008l : lMuzerolle et al.ll2005 



mm. 
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2. Methods 



We use a m odified version of the one-dimensional, two-zone disk model previously introduced in 



Zhu et all (|2010al lbl). including changes to the infall model, enhanced disk heating, and viscosity in the dead 
zone. Here we review our scheme. 



2.1. Surface density evolution 



The surface density of a disk is evolved based on the mass and angular momentum conservation equa- 
tions in cylindrical coordinates. 



dMi 



and 



d 



d 



d 



27rR^(m^n)-—iMiR^n) = 2TT—{R^WR^,i) + 27rAiiR,t), 



(1) 



(2) 



where S,- is the surface density, 17 is the angular frequency, M,- is the radial mass flux, Ws^ = RTjiVidQ/dR, 
and Vi is the viscosity. The subscript / denotes either the active layer ("a") or dead zone ("d"). The 
terms 2-Kgi{R^t) and 27r Ai{R,t) are the mass and an gular momentum flux per unit distance of infall material 
from an envelope cloud (|Cassen & Moosmanlll98lh . Then, assuming instantaneous centrifugal balance (see 



Zhu et al.ll2010br) . equations ([TJ and ^ can be simplified to 

Mi = 6^R'l^—{R'I^T.iUi)+ ; -Att — - {Ki{R,t)-gi{R,t)R^^{R)\ 

oR Mr at \GMrJ 



(3) 



where Mr is the sum of the mass of the central star and the disk mass within R. We then sequentially solve 
Equations Q and (d) to evolve the disk. 



In 



Zhu et al 



(1201 Obi) mass and angular momentum were added to the disk using the infall model of 
Cassen & MoosmanI (Il98lh . 

M. R^-'l^ 



g{R,t) = 



1 

R 



AttR, 
g{R,t) = OiiR>Rc, 



ifR<Rr 



and 



A{R,t) = g{R,t)R 



GM, 



1/2 



if R<Rr 



AiR,t) = if R>Rc 



(4) 
(5) 

(6) 
(V) 



Here Rc is the centrifugal radius, that is the outer radius at which mass is add ed to the disk at time t, and Mm = 
0.975c^/G is a constant total infall mass rate at a given cloud temperature ( Shull 19771 : iTerebey et al.lll984l) . 
However, this model has a singularity at R = R^; with finite grids this potentially causes non-convergent 
behavior at different grid resolutions. To avoid this we modified the infall model to eliminate the singularity 
by using a constant mass flux per unit distance. 
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Another problem occurs at the inner boundary, which is difficult to make very small because for nu- 
merical reasons (short time steps, dust evaporation, etc.) as well as fundamental uncertainties; for example, 
is the disk truncated by a stellar magnetosphere, and if so where, and is there an outflow from the inner 
disk edge. In addition, mass infall right at the inner boundary produces different results depending upon on 
precisely which inner boundary radius we choose. To minimize these problems we take inner boundary radii 
which are relatively small but still comfortably outside the expected point of magnetospheric truncation. We 
further assume that envelope material does not fall onto the disk inside 0.2Rc, justified on the basis that the 
well-known emergence of jets and outflows seen in even the earlies t protostellar phases sho uld prevent the 
lowest-angular momentum material from reaching the disk or star (IReipurth & Ballyll200lh . By tying the 
inner radius of infall to Rc we effectively assume that the same streamline denotes the boundary between 
outflow and inflow, such that the outflow cone retains the same opening angle, in this case, a half-angle of 
26.6°. This particular choice of opening angle is arbitrary and adopted mainly for numerical convenience. 



The mass infall rate of the modified model is 



if 0.2/?,. <R<Rc 



and 



g{R, t) = OifR< 0.2/?,. or /? > /?,.. 
The corresponding angular momentum added to the disk per unit distance is 



A(R,t) = g(R,t)R 



1/2 



if 0.2Rr <R<Rr 



and 



A(/?, = if /? < 0.2/?c or /? > /?,. 



(8) 



(9) 



(10) 



(11) 



A comparison of the radial mass infall profile of our model to that of the lCassen & MoosmanI (|l98ll ) model 
is presented in Figure E In total, q i ir mod el adds 1 1 % less angular momentum to the disk per unit mass 



infall than the 



increasing grid resolution. 



Cassen & MoosmanI (|l98lh model. These modifications result in better convergence with 



2.2. Temperature evolution 

The disk layer temperatures are determined by the balance between heating and radiative cooUng. For 
the active layer, the energy equation is 

C'£,,adtTa = 2heat,a ^ Gcool.a 

~ 2vis,a ~^ Ginfall.a ~^ 2grav,a ~^ ( ^ext~i ! ? 
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where Cs.a = T^acl^/Ta is the heat capacity of the active layer. Here Text characterizes the heating flux due 
to the irradiation of the disk by the stellar and (inner disk) accretion luminosity, and and are the optical 
depths of the active layer and the dead zone, respectively, 

Ta = ]^^aK{pa-,Ta) (13) 

and 

Td = -T,dKip d ,Td) (14) 

using the Rosseland mean opacity k taken from lZhu et al.l(l2009ah . In Equation (fT2l ). the first three terms aie 
local heating of the active layer due to the viscosity, the infall, and the gravitational potential energy change, 
respectively. The fourth term consists of the external heating (see below) and the radiative heating from the 
underlying dead zone. The last term includes the radiative cooling toward each side of the active layer. 

The energy equation of the dead zone is similar to that of the active layer, 

C-E^ddtTd = Qhea.t,d-Qcool,d^ (15) 

where Cs j = Tidcj^/Td is the heat capacity of the dead zone. If the active layer is optically thick, the energy 
equation is 



CT..ddtTd = Qvis.d + QmMi,d + Qgmv,d + -:r(^ ^ ^ - —aTf (16) 

3 1 + rj 3 1 + rj 



On the other hand, if the active layer is optically thin, the incident flux from the outside of the dead zone 
would be cr{TaT^ + T^^t) th^t energy equation is 



CT.,ddtTd = Qvis,d + Qial>in,d + Qgra.v,d + —(^{TaT^ + T^^t)—^-—(^T^-—^- (17) 

1+r^ 3 1 + T^, 



Again, the first three terms in equations ([T6] l and (iTTT i represent local heating, while the last two terms 
account for radiative heating from the outside of the dead zone and the radiative cooling. 

In the energy equations, the viscous heating is 

evis,- = ^WR^jn, (18) 

where Wr,^., = (3/2)E,f;0 and = a,c^,/rj. The viscosity parameter a,- is explained in detail in the next 
section. 

During infall the added material has smaller specific angular momentum than the disk material at the 
same radius. This results in a readjustment of the disk such a way that material moves inward. As we 
are assuming effectively instantaneous centrifugal balance, the increase in the gravitational potential energy 
driven by the readjustment process must be accompanied by the corresponding energy release. Here we 
assume that this heats the active layer only (2infaii,rf = 0), as this is the material directly impacted by the 
infalling matter. The heating by infalling material is then 



GMMin 3-2 JiR Rc) 

Ginfall a = ^ ' ' if 0.2/?, < R < (19) 
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and 



QinMla = Oif R<0.2Rc or R>Rc. 



(20) 



As accretion proceeds, the central stellar mass increases and the disk gravitational potential energy 
will become more negative. In response, even in the absence of viscosity disk material will move inward, 
implying additional accretion luminosity. The heating by this effect is 



Q 



grav,r ' 



2R 



(21) 



where M^, is change in the mass of the central star. 
The irradiation flux can be written as 



rp4 _ f*l^* ./acc^acc rp/\ 



(22) 



where L* and Lace are the stellar luminosity and the accretion luminosity, respectively, and Tenv is the en- 
velope cloud temperature. The coefficients and /ac e account f or the non-normal irradiation of the disk 
surface. For the stellar irradiation we use /* = 0.1 as in lZhu et all (l2010al) and assume the stellar luminosity 
follows the mass-luminosity relation 

which is an approximate power-law fit to pre-main sequence stars in the Taurus molecular cloud, using 
the luminosities and effective temperatures from Kenyon & Hartmann (1995), and adopting the Siess et al. 
(2000) eyqlutiona ry tracks to obtain the masses. The mass-luminosity relation is slightly modified from 
Zhuetal.1 (l2010bl\ 



log 



= 0.20+ 1.74 log — 



(23) 



The inclusion of external heating by inner disk accretion is another new feature of our calculations. 
The accretion luminosity is calculated as 

GM.M ^^^^ 



^acc — 



IR 







where we assume a typical T Tauri stellar radius. In this case the appropriate value of face is quite uncertain. 
At low to moderate accretion rates, magnetospheric accretion onto the star can occur at high latitudes, so 
that adoption of face = 0.1, similar to that used for the stellar photospheric irradiation, seems reasonable. On 
the other hand , at high accretion rates, the spectra of FU Ori objects provide no indication of magnetospheric 
accretion (Hartmann & Kenyon 19961) . and irradiation of the outer disk by a relatively flat inner disk should 
be much less effective (lBelllll999h . However, if a substantial infalling envelope surroun ds the disk, it can 
capture much of the accretion luminosity and reradiate a significant part toward the disk (lNattalll993h . We 
therefore use both /acc = 0.1 and 0.01 to examine the importance of this heating. 



The last term in Equation (1221 ) is the flux from the envelope cloud whose temperature is assumed to 
20 K. Thus, the stellar luminosity irradiation Q^^ and the accretion luminosity irradiation Qi^^ on the active 
layer become 

16 f^L^ Td 



Q* = 



3 47r/?2 l + ri 



(25) 
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and 

el6 /acc^acc Ta t'^c\ 
''-^y^^TT^- ^^^^ 

We note that the accretion luminosity irradiation 2acc should be distinguished from the local viscous ac- 
cretion heating 2vis- The relative importance of the individual heating terms, and their effects, during disk 
evolution will be discussed in 01 



2.3. Disk Viscosity 



The viscosity parameter a, is the sum of the MRI viscosity parameter a^,/ and the GI viscosity pa- 
rameter aqj. The MRI viscosity parameter is assumed to have a fixed value of aum only if a region can 
sustain the MRI. Thus, the active layer viscosity parameter is always set to aM,a - olmri while the dead zone 
has MRI viscosity only if the midplane temperature is higher than a critical temperature Tmri to produce 
sufficient ionization levels. For the dead zone, we consider a residual viscosity as well as MRI viscosity and 
GI viscosity, = aM,rf + 02,£/ + ard- 

The idea of the dead zone residual viscosity (DZRV) ard is based on recent numerical magnetohy- 
drodynamic simulations suggesting that magnetic turbulence in the active layers can drive h ydrodynamic 
turbu l ence in the dead zone, irnplying a non-zero r esidual viscosity parameter ~ 10~^ - 10"^ (IBai & Stone 



201ll : lQkuzumi & Hirosell2011 



Gressel et al.ll2012h . Thus, for non-zero DZRV model we set 



Q!rd = min 10 , /rdOMi?/ 



(27) 



where /rd is the efficiency of accretion in the dead zone whose value is chosen to be < 1 ; this is intended 
to limit the effect of the active-layer induced turbulence such that the mass accretion rate of the dead zone 
(approximately) does not exceed that of the active layer (M^ < M^). This seems intuitively reasonable. We 
consider the upper limit f^d = 1 and consider a case with f^d = 0. 1 as it is unlikely that the active layer can be 
that effective in driving accretion. 



Finally, the GI viscosity parameter is the same as in lZhu et al.l (l2010ah . 



(28) 



where Q is the Toomre parameter. 



3. Results 
3.1. Initial conditions 

We start with a 0.1 M© central protostar surrounded by an Mc = 1 M© cloud. We parameterize the 
cloud rotation in terms of a; = Qc/^b, where Qc is the (constant) angular frequency of the initial cloud, and 
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= is the breakup angular frequency at the outer cloud edge, and c.v is the (uniform) cloud 

sound speed. Our fiducial models assu me to = 0.03, whic h results in 15 % larger cloud angular frequency 
than that used in the fiducial model of IZhu et all (1201 Obi) ($7^ ~ 1.15 x 10~^^ rad s"^ in our model). We set 
the maximum non-thermally ionized surface density to 100 g cm"^ for our fiducial choice and assume 
it is constant. We assume Tmri = 1500 K and omri = 01 for all calculations. We adopt a cloud envelope 
temperature of T^nv = 20 K, which yields a constant infall rate of ~ 3.4 X lO-^Movr-^ . This is 20 % smaller 
than the infall rate for conventional singular isothermal collapse model (|Shull 19770 because of our modified 
infall model (see 32.11 ). The infall lasts for ~ 0.24 Myr, adding 0.8 M© to the central star + disk in total. 



3.2. Zero dead zone viscosity model 



Zhu et al.l (|2010al ) found that GI moves matter from the outer disk to the inner disk, leading to a pileup 



at /? ~ 2 AU because GI is increasingly ineffective at small radii. Eventually enough material piles up to trap 
thermal energy that m akes > Tmri, tu rns on the MRI thermally in the dead zone, and thus produces an 
outburst of accretion. IZhu et al.l (|2010b|) investigated the long-term evolution of such disks and found that 
the evolution can be divided into three stages. The evolution starts with a quasi-steady disk accretion since 
the infall is to small radii where the inner disk can become hot enough to sustain the MRI thermally. Then, it 
turns into the outburst stage as the infall occurs at radii > 1 AU. After infall stops, the disk enters the T Tauri 
phase, having only a few Gl-driven outside-in outbursts with a low mass accretion rate in between bursts. 

Figure El shows the mass accretion rate and the mass of the central star, the disk, and the central star -i- 
disk as a function of time. The modified infall model and additional heat ing sources discussed in §2 produce 
no qualitative difference in the overall evolution from Izhu et al. (2010b). 



The new heating sources for the first 0.3 Myr of the evolution at R = I and 10 AU are presented in 
Figures |3]and|4j together with the previously considered sour ces. In these figure s, newly added terms in this 
paper are plotted in color while other terms have considered in lZhu et all (|2010ah and drawn in black. As one 
can see, heating by the change in gravitational energy is usually several orders of magnitude smaller than 
other terms so that it makes no change in the evolution. Infall heating provides comparable amount of heat 
to the active layer but is only limited to the region material falls onto, increasing the local disk temperature 
slightly. However, the large increase in the accretion luminosity during outburst produces enough irradiation 
to make the outer disk temperature increase dramatically. This is shown in Figure|5] where we show the mass 
accretion rate during a single outburst and the radial profiles of the disk surface density and the midplane 
temperature before outburst, at the maximum accretion rate, and at the end of the outburst. The temperature 
increase at the outer disk during outbursts does not affect the long-term evolution, because the viscous time 
of the outer regions (~ 10^ yr) is much longer than the outburst timescale of ~ 10^ yr. 

We also call attention to the jump in temperature at radii < 0.5 AU, which is due to thermal instability 
dZhu et al.ll2010ar) . This increase in temperature, which affects the behavior of the outburst of accreting 
material onto the central star, would not have been found if we had taken an inner radius of > 1 AU (see 
discussion in §4). 



-9- 



The Gl-driven outside-in bursts accrete ~ 0.027 Mq of material onto the central star and last for ~ 1340 
years on average. Although the duration obtained in our calculation is longer than the typical outburst 
timescale seen in F U Qri, the timesc ale of an outburst can be scaled with a choice of 

'^MRI since ^ 

R^/uoc ttMRi"^ (see lzhu et al. IboiOal ). The outside-in bursts are rare, because the disk needs a lot of material 
to trigger MRI through GI while it is difficult to do so with zero DZRV. The duty cycle of this model is thus 
pretty small, ~ 0.06 during outburst stage and ~ 0.005 during T Tauri phase. The details of disk properties 
at 0.24 and 1 Myr and of outside-in bursts are summarized in Table[T] All the outburst quantities in the table 
are time-averaged values after the initial quasi-accretion phase while those vary with time. 



3.3. Non-zero dead zone residual viscosity model 

While the modified infall model and additional heating sources with zero DZRV make no qualitative 
change in the overall evolution, a finite DZRV makes a lot of difference in the long-term evolution and in 
the single outburst behavior as well. Figure |6]shows the mass accretion rate and the mass of the central star, 
the disk, and the central star -i- disk as a function of time. During infall the disk has a quasi-steady accretion 
phase at the beginning (t < 0.08 Myr) and the outburst stage follows, as in the zero DZRV case. However, 
the evolution after the accretion phase is different in that the non-zero DZRV model shows a lot of smaller 
outbursts instead of a few large outbursts. 

Figure [T] and [8] show the heating sources of active layer and dead zone of the non-zero DZRV model 
at /? = 1 and 10 AU, respectively. The major difference between this model and the zero DZRV model is 
that dead zone viscous heating provides a significant amount of heat at the inner disk even after infall ends. 
We have run calculations with different Tmri (1300 K and 1800 K) and found that the overall features are 
not sensitive to the choice of Tmri. At the outer disk, the accretion luminosity irradiation is still important 
during bursts while the temperature increase is not as dramatic as in the zero DZRV model due to lower 
accretion peak. Infall and gravitational heating are less important than others. 

Figure |9t a) shows the mass accretion rate during a single outburst with non-zero DZRV. Radial surface 
density and midplane temperature profiles at the beginning, at the maximum accretion rate, and at the end of 
the burst are presented in Figure 0b) and (c). The accretion behavior during a single outburst is remarkably 
different from that of the zero DZRV model. The outburst has a peak of Mmax ~ 10"^ MQyr~', which is 
about two orders of magnitude smaller than that of the outside-in bursts. In addition, the accretion rate 
initially shows a rapid increase but has a slow rise time to its peak and a slow decrease after the peak as 
well. In this model, the dead zone is able to transport material with the help of the non-zero DZRV. Thus, 
the inner disk can be heated viscously and outbursts are initiated at the inner boundary of the disk before 
the GI piles up enough material at the middle of the disk (/? ~ 2 AU) to initiate the MRI, which is the case 
for outside-in bursts. The ionization front propagates out to several AU from the inner boundary. Since 
the inside-out bursts have an MRI active inner boundary from their initiation, the disk continues to dump 
material from its innermost part during the whole bursts. Therefore, the system is not able to show a huge 
accretion rate as seen in outside-in bursts, but only generates moderate accretion rate. Note again that the 
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outburst triggers first at small radii, inside of ~ 0.5 AU (§4). 

On average, the mass accreted onto the central star during a single inside-out burst is ~ 1.5 x 10"^ Mq 
and it lasts ~ 450 years. The inside-out bursts occur frequently enough to get a duty cycle of ~ 0.16 during 
outburst stage and ~ 0.06 during T Tauri phase. The disk properties and outburst details of the non-zero 
DZRV model are summarized in Table [T] 



3.4. Efficiency of accretion luminosity irradiation 

As shown in the previous sections, irradiation by the inner disk plays an important role during outburst 
on the temperature profile at the outer disk. However, the efficiency of the accretion luminosity irradiation 
is uncertain as far as the non-normal irradiation of the disk is considered. We thus test the effect of changing 
/acc to 0.01, which is ten times smaller than the fiducial value. 

We found essentially no change in the overall evolution of both zero and non-zero DZRV cases, since 
the accretion luminosity irradiation is several orders of magnitude smaller than main heating sources - active 
layer viscous heating and stellar irradiation - during the quiescent phase. During outbursts, however, the 
accretion luminosity irradiation still dominates the heating even with a ten times smaller efficiency, making 
a significant difference to the outer disk temperature. Not surprisingly, the increase in outer disk temperature 
during bursts is smaller than the standard cases by a factor of ~ 2. 



3.5. Accretion Efficiency in Dead Zone 

Intuitively, it seems unlikely that the turbule nce generated by the ac tive layers within the dead zone can 



transport as much mass as the active layer (e.g., iHartmann et al.ll2006h . In our models this happens when 
"rd ^ 10~^, where > 10^ g cm~^. This could be an overestimate of the efficiency with which the MRI 
turbulence in the active layers drives accretion in the dead zone. We therefore adjust the dead zone accretion 
efficiency /rd to 0.1 so that dead zone only has an accretion rate of ~ 10 % of the active layer at most. 

Figure [To] shows the mass accretion rate and the mass of the central star, the disk, and the central star 
-I- disk as a function of time. Initially, the evolution resembles that of the standard zero DZRV model rather 
than the non-zero DZRV model; the system shows a distinct outburst phase during infall. This is because 
the mass that the dead zone can carry is now limited and thus generates less viscous heating at small radii 
than the standard non-zero DZRV case. Therefore, mass piles up at large radii through the GI before inner 
disk gets heated and triggers inside-out bursts. After infall ends, however, we still see inside-out bursts with 
much less frequency than the standard non-zero DZRV model, which is again due to less viscous heating 
at the inner disk. The duty cycle during T Tauri phase of this model is only 0.015, which is four times 
smaller than that of the standard non-zero DZRV model. First three outbursts after infall ends are outside-in 
bursts, since disk already collects enough material at outer disk during infall to make them. This emphasizes 
importance of understanding the effect of MRI turbulence on dead zones (§4). 
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3.6. Dependence on Ea and umri 



While we use = 100 g cm ^ as our fiducia l value, sev eral stud ies have pointed out that the active 



layer more likely has a lower surface density (e.g. ISano et al. l 2000; Ba i & Goodmanll2009l) . We thus test 
Sa = 20 g cm"^ and adjust the MRI viscosity parameter umri = 0.05 to maintain roughly the same mass 
accretion rate (M oc aS) during the quiescent phase as the standard cases (and also in agreement with 
typical T Tauri accretion rates). 

Figure [TT] shows the mass accretion rates of both zero and non-zero DZRV models adopting the lower 
value of Sa. The mass accretion rate during a single outburst and radial surface density and midplane 
temperature profiles at the beginning, at the maximum accretion rate, and at the end of a single outburst 
of the both models are presented in Figure [12] Since the outburst timescale depends on the MRI viscosity 
parameter (Afburst oc aMR/~^)> the details of the outbursts, such as outburst duration and peak accretion rate, 
vary. However, the overall evolution as well as the initiation of outbursts remain the same. We see GI- 
induced MRI-driven outside-in bursts in the zero DZRV case and viscously triggered inside-out bursts in the 
non-zero DZRV case. This is because the overall evolution and the initiation of outbursts are governed by 
the mass accretion during the quiescent phase, whi ch we manage t o be unchanged. We note that the shorter 
timescales are in better agreement with FU Ori (see lZhu et al The disk properties and outburst details 

are summarized in Table [I] 



4. Discussion 



Our simulations show that it is possible to obtain inside-out triggering of accretion outbursts as well 
as outside-in bursts (e.g., Zhu et al. 2010a,b), with the former enhanced if there is finite dead zone residual 
viscosity. The two types of outbursts were also obtained in the model developed by Bell & Lin (1994; BL) 
for FU Ori outbursts. In the BL model, the outbursts were due to thermal instability (TI), plus an assumed 
increase in a from a very low value to a much higher value. As Zhu et al. (2007, 2008) showed, the TI 
model is inconsistent with observations of FU Ori, because the high temperatures required limit the region 
of rapid accretion to smaller radii than inferred from modeling the spectral energy distribution including 
Spitzer Space Telescope data. Nevertheless, the finite dead zone residual viscosity models are qualitatively 
similar to the basic feature of the BL models which produce inside-out bursts; a small but finite viscosity 
allows material in the inner disk to produce enough trapping of viscously-generated heat to trigger a higher 
viscosity and eventually an outburst. As BL showed, such inner disk triggering leads to outbursts with slow 
rise times, qualitatively consistent with the observed outburst of V1515 Cyg (Herbig 1977; M. Ibrahimov, 
personal communication). 



BL showed that outbursts with rapid rise times, such as observed in FU Ori and V1057 Cyg (IHerbig 



19771) . required outside-in accretion events. In the BL model, a large outer perturbation of the disk was 
required. In modem models, the event is triggered by GI, which piles up material at larger radii than po s 



sible in the TI mod el (|Armitage et al.ll2001 



Vorobvov & Basu 



Martin et al.ll2012alJb|) . Our current results build upon those of 



2006, 



2007 



Zhu et al 



2009L 120101 : IZhu et al.ll2010al lbl: 
m that we clearly identify 
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some inside-out bursts during the main phase of infall (they were actually present in the Zhu et al. simulation 
as well but were not emphasized). 

Our results also bear similarities to outbursts in models of cataclysmic variables (CV). Two different 
types of outb ursts (outside-in and inside-out) in accretion disks around dwarf novae were first predicted by 
SmakI (|l984l) . In CV models, mass transfer from the secondary rises the effective temperature of a disk 
annulus to 5000 - 8000 K, which corresponds to hydrogen recombination inside the disk and thus triggers 
the TI (lMenouetalJll999h . While the dead zone residual viscosity parameter is the feature that changes 
outburst behavior in our model, CV models use the mass transfer rate to generate two different outbursts. If 
the mass transfer rate is high, the accumulation timescale of transferred material is shorter than the viscous 
timescale, allowing material piles up at outer disk. Thus, outside-in bursts are triggered. In contrast, if 
mass transfer rate is low, the accumulation timescale becomes longer than the viscous timescale so that 
the outbursts outside-in outbursts are replaced by inside-out ones. The resulting outburst behaviors of CV 
models are similar to ours; inside-out bursts have a small er accretion peak and a slower rise time than 
outside-in ones (see Figures 3 and 4 of iHameury et al.lll998[) . 



It is worth emphasizing that the outburst behavior of systems, observed at optical and near-infrared 
wavelengths, is a result of accretion onto or near the star, i.e. at radial scales < 0.1 AU. While our models 
do not reach magnetospheric or stellar radii, our inner boundary radius of 0.2 AU is small enough to capture 
behavior (thermal instability, inside-out bursts) which cannot be seen in simulations with inner boundaries 
> 1 AU. Thus, while our own treatment of non-steady accretion has its limitations, time histories of accretion 
in simulations with large inner disk r adii must be treated with special cauti on. Using an inner boundary of a 
few AU, as in the series of papers by lvorobyov & Basul (120061 boOvL 120091) or bunham & Vorobvoyl(l2012l) . 
one would not find t he outburst behavior characterized by our models or those of Armitage et al. ( 200ll) or 
Martin et al. (2012a b). 



Along these lines, we have found that the precise duration of the inside-out bursts during infall is 
sensitive function of the value of the inner radius. We do not explore this further here because there are 
other major uncertainties in our treatment, such as the assumption of a constant active layer surface density, 
the characterization of the GI via an a viscosity, and the way in which we implement a finite dead zone 
viscosity. In the following paragraphs we discuss these issues in turn. 

Advanced MHD treatments of the active layer are complex and involve a number of unknowns, such 
as whether low-energy cosmic rays can penetrate the accretion-driven winds, grain growth and settling, 
the presence of metal ions, etc. (e.g., Sano et al. 2000; Ilgner & Nelson 2006, 2008; Hirose & Turner 
2011; Perez-Becker & Chiang 2011). Martin et al. (2012a,b) argue that a large critical Reynolds number 
is necessary for transport to occur, resulting "active" regions very different than the constant we use. 
However, Martin et al. (2012b) predict essentially no accretion onto the central star in between outbursts, 
whereas the pre-outburst spectrum of V1057 Cyg shows emission lines typical of T Tauri stars accreting at 
~ lO~^M0yr"^ (Herbig 1977). More recently. Miller et al. (2011) showed that the classical (accreting) T 
Tauri star LkHa 188-G4 underwent an FU Ori-type eruption in 2009. Finally, a pre-outburst spectrum of the 
FU Ori object V733 Cep (Reipurth et al. 2007) also showed characteristic accreting T Tauri emission lines 
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(B. Reipurth, personal communication). Thus the pre-outburst state of at least some FU Ori objects is one 
of accretion at rates typical of T Tauri stars, which we obtain with our adopted value of Q!qSq, at least in the 
inner disk. 

Similarly, our treatment of the GI with an a viscosity is crude; as in the case of the MRI, three- 
dimensional simulations are required to treat the GI properly (e.g.. Rice et al. 2003; Boley et al. 2006; 
Durisen et al. 2007, and references therein). Two-dimensional simulations also do a better job of capturing 
the GI than our treatment (e.g., Vorobyov & Basu 2006, 2007, 2009, 2012). However, as pointed out by Zhu 
et al. (2009), and as shown in Zhu et al. (2010a,b) and the present simulations, the GI becomes harder and 
harder to sustain as one moves to smaller radii, whereas triggering of the MRI becomes easier. 

Finally, the presence and behavior of non-zero viscosity in dead zones is highly uncertain. While 
simulations of resistive, stratified disks in shearing boxes appear to show that the active layer produces a 
non-zero a in the central "dead" layers (Fleming & Stone 2003; Okuzumi & Hirose 2011; Gressel, Nelson, 
& Turner 2012), the level to which this occurs, the amount of mass transport involved, and precisely where 
energy is dissipated is unclear. For example, our model assumes that the wave energy is dissipated near the 
midplane, and thus the heat generated can be trapped radiatively by the opacity of the disk at the midplane; 
however, it is possible that dissipation is concentrated at higher levels (N. Turner, personal communication). 

In all, these uncertainties show that our current results must be taken as suggestive rather than predictive 
for variations of accretion in young stellar objects. Nevertheless, these simple models, which can be evolved 
easily for significant evolutionary timescales, illustrate the potential information on transport processes in 
protostellar and protoplanetary disks that might ultimately be gleaned from the observed accretion outbursts. 
Due to the increasing monitoring of young stellar objects, it is becoming increasingly clear that a wide 
variety of accretion behavior is exhibited in young stars, emphasizing that further progress on challenging 
problems of globally simulating MRI and GI in protostellar disks may pay rich dividends. 



5. Summary 



In this paper, we have extended the one-dimensi onal, two-zone mo del of long-term protostellar disk 
evolution with infall, which is previously introduced in lZhu et al.l fcOlOal ibh. Our modified models include a 



revised treatment of infall, enhanced disk heating, and possible non-zero viscosity in th e dead zone. Whil e 
the former two changes produce no qualitative difference in the overall evolution from the lZhu et al.l(l2010bh . 
we find that the presence of a small but finite dead zone viscosity can trigger inside-out bursts initiated at 
or near the inner edge of the disk through dead zone viscous heating, instead of Gl-induced MRI-driven 
outside-in bursts with zero dead zone viscosity. These inside-out bursts not only bear a qualitative resem- 
blance to the outburst behavior of one FU Ori objects, VI515 Cyg, but emphasize a careful treatment of the 
inner disk regions in simulations. 

Given the uncertainties, our results are rather suggestive than predictive. However, two types of out- 
bursts seen in FU Ori objects can be successfully reproduced by the simple a treatment in the dead zone. 
This difference in accretion behavior could be a potential probe of transport efficiency in the dead zone. 
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Fig. 1. — The mass infall rate at ? = 0.05 Myr of ICassen & MoosmanI (119811) mode l (asterisks) and our 
consta nt g{R) model (diamonds). The centrifugal radius at the time is Rc ~ 1.24 AU. ICassen & Moosman 
(1198 ih model gives ~ 1.3 x 10"^ Mq yr~' of mass infall rate at the nearest grid to Rc, which is about an 
order of magnitude large to be fitted in the figure. The mass infall rate inside 0.2 Rc is set to zero in our 
model to imitate protostellar outflows (see text). Dotted curves show analytic estimates of the mass infall 
rate of each model. 
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Fig. 2. — (a) Mass accretion rate and (b) mass of the central star + disk (dotted curve), mass of the central 
star (solid curve), and mass of the disk (dashed curve) w^ith time for the standard zero DZRV model. 
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Fig. 3. — Various heating sources of the active layer (upper panels) and the dead zone (lower panels) at 
R = I AU for the tirst 0.3 Myr (left panels) and during a single outburst (right panels) with zero DZRV. 
Newly added hea ting sources in this paper are plotted in color while the heating sources have considered in 
dzhu et allboiOal lbl) are presented with black curves. Upper panels: 2vis,a (black solid), 2* (black dotted), 
Ginfaii (red dashed), 2acc(green dashed), and Qgrav,a (blue dashed) are presented. Lower panels: 2vis,t/ (black 
solid) and 2grav,d (blue dashed) ai-e presented. 
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Fig. 4.— Same as Figure|3]but at /? = 10 AU. 
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Fig. 5. — (a) Mass accretion rate of a single outburst of the standard zero DZRV model. The " drop out" 



in acc retion during the middle of the outburst is an artifact of the one-dimensional treatment (|Zhu et al 



2010ah . Radial profiles of (b) the surface densities and (c) the midplane temperatures before the outburst 



(solid curves), at the maximum accretion rate (dashed curves), and at the end (dash-dotted curves) of the 
burst are plotted. In panel (b), the curves with higher surface densities at the inner disk represent the dead 
zone surface density while the lower ones extend further out represent the active layer. In panel (c), the 
dotted horizontal line represents the MRI activation temperature Tmri = 1500 K. 




Fig. 6. — (a) Mass accretion rate and (b) mass of the central star + disk (dotted curve), mass of the central 
star (solid curve), and mass of the disk (dashed curve) with time for the standard non-zero DZRV model. 
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Fig. 7. — Various heating sources of the active layer (upper panels) and the dead zone (lower panels) at 
/? = 1 AU for the first 0.3 Myr (left panels) and during a single outburst (right panels) with non-zero DZRV. 
Newly added hea ting sources in this paper are plotted in color while the heating sources have considered in 
dzhu et allboiOal lbl) are presented with black curves. Upper panels: 2vis,a (black solid), 2* (black dotted), 
Ginfaii (red dashed), 2acc(green dashed), and Qgrav,a (blue dashed) are presented. Lower panels: 2vis,t/ (black 
solid) and 2grav,d (blue dashed) ai-e presented. 
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Fig. 8.— Same as Figure|7]but at /? = 10 AU. 
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Fig. 9. — (a) Mass accretion rate of a single outburst of the standard non-zero DZRV model. Radial profiles 
of (b) the surface densities and (c) the midplane temperatures at the beginning (solid curves), at the maximum 
accretion rate (dashed curves), and at the end (dash-dotted curves) of the burst. In panel (b), the curves with 
higher surface densities at the inner disk represent the dead zone surface density while the lower ones extend 
further out represent the active layer. In panel (c), the dotted horizontal Une represents the MRI activation 
temperature Tmri = 1500 K. 




Fig. 10. — (a) Mass accretion rate and (b) mass of the central star + disk (dotted curve), mass of the central 
star (solid curve), and mass of the disk (dashed curve) of non-zero DZRV model as a function of time, with 
10 % of accretion efficiency /rd in the dead zone. 



cj> -10 



© 



cn 
O 



0.0 



0.0 



0.5 



1 .0 



1 .5 



time [Myr] 




0.5 



1 .0 



1 .5 



time [Myr] 



Fig. 1 1 . — Mass accretion rate of zero DZRV model (upper) and non-zero DZRV model (lower) as a function 
of time, with Ea = 20 g cm~^ and aMRi = 0.05. 
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Table 1 . Parameters and results 





"MRI 












Mm" 




/M'max" 




^/burst'' 


Do' 








(g cm ^) 


(AU) 


(Mq) 


(Mq) 


(Mq) 


(Mq) 


(Af©) 


(Mq yr-1 


) 


(yr) 






zero 


0.01 


100 


29/5.5 


0.64/0.76 


0.26/0.14 


0.21/0.08 


0.05/0.06 


0.027 


1.11 X 10" 


-3 


1340 


0.06 


0.005 


non-zero 


0.01 


100 


24/6.5 


0.65/0.80 


0.25/0.10 


0.20/0.04 


0.05/0.06 


1.5 X 10"^ 


3.17 X 10" 


-6 


450 


0.16 


0.06 


zero 


0.05 


20 


33.9/7.6 


0.65/0.80 


0.25/0.10 


0.22/0.08 


0.03/0.02 


0.028 


4.34 X 10" 


-3 


310 


0.016 


0.002 


non-zero 


0.05 


20 


33.9/7.6 


0.65/0.82 


0.25/0.08 


0.22/0.05 


0.04/0.03 


4.6 X 10"^ 


4.26 X 10" 


-5 


200 


0.08 


0.015 



^Quantities are taken at 0.24 and 1 Myr. 

''Outburst quantities are averaged over time after infall ends. 

'^Duty cycle during outburst stage. 
*^Duty cycle during T Tauri phase. 



